Tracking heterogeneous dynamics during the ^-relaxation of a simple glass-former 
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(Dated:) 

We study the relaxation process in a simple glass-former - the KA lattice gas model. We show that, 
for this model, structural relaxation is due to slow percolation of regions of co-operatively moving 
particles, which leads to heterogeneous dynamics of the system. We find that the size distribution of 
these regions is given by a power-law and that their formation is encoded in the initial structure of 
the particles, with the memory of initial configuration increasingly retained with increasing density. 

PACS numbers: 



Glass-forming systems are disordered materials whose 
relaxation dynamics becomes extremely slow on decreas- 
ing the temperature or increasing the density. Despite 
intensive research in this domain, a proper understand- 
ing of these materials is still missing [l|, |2| . The observa- 
tion, in experiments and simulations, of spatio-temporal 
dynamical heterogeneity in glass-formers has been an im- 
portant step forward in elucidating the mechanism for re- 
laxation processes and currently considerable research is 
being focused on understanding the lifetime and spatial 
extent of these dynamical heterogeneities. 

A key question in this domain is the conundrum re- 
garding whether or not there is a causal link between 
structural properties of glass-forming systems and their 
dynamical behavior. It has recently been shown Q that 
the spatial heterogeneity in the propensity of particles 
to move is correlated with the local environment of the 
particles, characterized by the local Debye- Waller fac- 
tor. However, this influence seems to exist only over 
timescales which are much less than the structural relax- 
ation time [J] . It was also observed that structural prop- 
erties are indeed correlated with collective dynamical 
fluctuations, but no quantitative analysis was made 
More recently, it has also been reported p that for a su- 
percooled liquid configuration, its localized low-frequency 
normal modes correlate with the irreversible structural 
reorganization of its constituent particles. 

In recent times, extensive studies of kinetically con- 
strained models (KCM-s) 0] , which are one of the sim- 
plest models showing glassy dynamics, have been car- 
ried out in order to understand their relaxation process. 
These models, which are motivated by the hypothesis 
that the slow dynamics in glass-formers is only due to ge- 
ometrical constraints, show heterogeneous dynamics sim- 
ilar to real glass-formers In this Letter, we stud 
one such KCM - the Kob-Andersen (KA) lattice gas 
in which particles are allowed to move on a lattice fol- 
lowing certain dynamical rules and which at high densi- 
ties shows si gna tures of apparently diverging relaxation 
times 0, Hoi I 111) . Recently it has been proven analyti- 
cally [IH and numerically [3] that for this model there 
exists no dynamical transition at finite density, p, and it 



was argued [12j that eventually, due to the presence of 
migrating macro- vacancies, the system relaxes, albeit ex- 
tremely slowly. However, from a practical point of view 
this model still is a good model for a glass-forming sys- 
tems. Using Monte Carlo simulations we show that its 
structural relaxation is related to the growth of mobile 
regions and that this process quickly slows down with in- 
creasing p, resulting in the observed slow dynamics. We 
also demonstrate that the formation of the mobile regions 
is directly related to structural properties of the system. 

We have studied the 3d version of the KA model: N 
particles populate a cubic lattice of size L 3 with the con- 
straint that a lattice site can be occupied by only one 
particle. All possible configurations have the same en- 
ergy and thus the same Boltzmann weight. The imposed 
stochastic dynamics consists of the following process: A 
randomly selected particle can move to any one of the 
neighboring empty lattice site provided it has m or fewer 
occupied nearest neighbor site and that the target empty 
site has to + 1 or fewer occupied nearest neighbor sites. 
A choice of m = 3 results in glassy dynamics for this 
model For efficient sampling of the configuration 
space at high a = N/L 3 , we have carried out event-driven 
Monte Carlo [9] simulations of the model. Using periodic 
boundary conditions, we have investigated system sizes 
L = 20, 30, and 50, which avoid finite size effects, with 
densities spanning from p = 0.65 to p = 0.89. 

Experiments and simulations in which the motion of 
single particles were tracked have helped to demonstrate 
the existence of heterogeneous dynamics in glassy sys- 
tems 0). E.g., by measuring the self part of the van Hove 
function G s (r 7 1), i.e. the distribution of particle displace- 
ments (G s (r,t) = (S(r — |r,-(t) — r,-(0)|)), where rj(t) de- 
notes the position of particle i at time t) , it has been pos- 
sible to demonstrate that the particles have varying mo- 
bilities. In Fig. [I] we show G s (r,t) for different times t, 
measured in units of Monte Carlo steps, at p = 0.87. For 
diffusive motion, G s (r,t) is a Gaussian and we observe 
that for the KA model, like other glass-formers this 
Gaussian behavior is only observed at times [t rj 5 x 10 8 ) 
that are much larger than the structural relaxation time 
r a (defined as the time at which the self-intermediate 



FIG. 1: G a {r,t) for different times at p = 0.87. The 
a— relaxation time at this density is r a ~ 1.8 x 10 7 . 

scattering function has decayed to 1/e, which at p = 0.87 
is r a k 1.8 x 10 7 ). At intermediates times, we see that 
G s (r,t) has an exponential tail, a signature of the pres- 
ence of rare events in the dynamics of the particles, sim- 
ilar to other glassy systems [l5j | . Thus we can conclude 
that while most of the particles remain frozen at their ini- 
tial positions (resulting in large values for G(0,t)), there 
is a small population which is extremely mobile, i.e. the 
system has a very heterogeneous dynamics. 

By tracking the mobile particles, i.e. particles which 
contribute to the tail of G s (r, t) , we see that at short 
times they explore a compact region (a "blob") around 
their initial locations. With increasing t these blobs 
slowly expand and coalesce with other blobs to form a 
labyrinthine structure. This structure allows particles, 
which were hitherto confined to one blob, to travel longer 
distances. In Fig. [5] we have plotted the lattice sites 
(marked by green and black spheres) visited by two such 
mobile particles at p = 0.88 and t — 10 7 (at this density, 
r a « 4.4 x 10 8 ). It can be clearly seen that the dynamics 
is spatially heterogeneous since the trajectories consist 
of connected blobs. Note that initially each blob is a re- 
gion of cooperative motion since all the particles which 
occupy these lattice points are found to be mobile. As 
the backbone is formed, particles from one blob begin to 
explore other blobs. The relaxation of the system thus 
happens with the slow growth of this backbone, within 
which the particles can move relatively quickly. This is 
demonstrated in Fig.O: indeed, the two trajectories over- 
lap even though they originate from two different lattice 
sites which initially did not belong to the same blob. 

Further insight into the spatial nature of the relaxation 
process can be obtained by observing the so-called "mo- 
bility regions" 

0, 

Il7l Il8l: A lattice site is defined to 



be an "active site" if either a particle or a vacancy has 
moved out of it during the time of observation and the 
collection of these sites constitute the mobility regions. 
Earlier studies of KCM-s have shown that the active sites 
tend to cluster and act as seeds for subsequent mobility 
[HI, 12 > EH and that it is possible to extract from the 



FIG. 2: Green and black points: Sites that have been visited 
by the trajectory of two mobile particles. Red blobs: The 
mobility regions, t — 10 7 and p — 0.88. 



mobility regions a lengthscale which increases with den- 
sity In Fig. M 

we have also included the location 
of the active sites and we can see that these lattice sites 
are indeed clustered in space (marked by the red blobby 
shapes) and have a labyrinthine structure. It gives us 
an idea of the structure of pathways, at high densities, 
available to the mobile particles for exploration. We can 
clearly see that, on this time-scale, the two mobile parti- 
cles have only explored a part of the available volume and 
that the geometry of the mobility regions and the blob- 
structure of the trajectories are intimately connected to 
each other. 

The number density of active sites, n act (t), allows us 
to estimate the volume accessible to the mobile parti- 
cles and in Fig. [3] we plot n act (t) for different p. We 
see that, at short times, n act (t) increases quickly and we 
find n ac t(t) ~ 1 — exp(— 1/9) with 9 w 5, independent 
of p. This regime corresponds to the initial growth of 
the blobs. Subsequently the shape of n act (t) depends 
strongly on p. For p — 0.80 the number of inactive sites, 
1 — n act , decays with a stretched exponential tail, with 
a stretching exponent of around 0.6, a functional form 
that is found for all p. At even larger p, n act (t) shows 
three regimes, with the second regime being a period of 
extremely slow growth, almost logarithmic and thus sim- 
ilar to the coarsening process in disordered media [l9j |. 
Note that at short t the typical distance between the 
blobs increases with p and the growth of the blobs slows 
down with increasing p, since it needs the presence of 
active sites (which are rare at high density). This is the 
reason why the increase of n act (t) at intermediate times 
becomes very slow with increasing p. We also observe 
that at p — 0.89, which is higher than the density of 
p c = 0.881 at which an apparent divergence of relaxation 
timescales was observed [9J, n act (t) is still an increasing 
function, suggesting that the system will eventually re- 
lax 12j. Note that although these three regimes are in 
qualitative agreement with the predictions of Ref. 12j |. 
there are important differences since, e.g., the diffusion 
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FIG. 3: Growth of number of active sites, function 
of time, for different p. The arrows mark r a . 

of macro- vacancies discussed in Ref. |l2[ would lead to a 
linear growth of n ac t(t) at long times, a behaviour which 
is not seen in Fig. [3] This might be due to the fact that 
the calculations presented in Ref. [13] apply only at den- 
sities that are extremely close to 1.0. 

In order to characterize the geometry of the growing 
clusters of active sites we have calculated P(s), the dis- 
tribution of clusters that have exactly size s. In Fig. 2J 
we have plotted P(s) for p = 0.88 and different times. 
At short times, t = 10, P(s) has an exponential shape. 
This corresponds to the initial geometry of active sites at 
few random locations when particles explore their neigh- 
borhood. With increasing time P(s) quickly transforms 
into a power-law, P(s) ~ s _ly , with an exponent v ks 1.6. 
This indicates that the growth process is different from 
random percolation for which v — 2.2. This difference is 
likely related to the fact that there is a wide variation in 
the size of the cooperatively rearranging regions seen at 
short times. At even later times, the largest cluster starts 
growing and the tail in P(s) shifts to larger and larger 
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FIG. 4: Distribution of cluster size, P(s), where s is the size 
of cluster, for different times at p = 0.88. Also shown (dashed 
line) is an exponential fit to P(s) at t — 10 and (solid line) 
the function s -1,6 to highlight the emergence of power-law 
behavior of P(s) at intermediate times. 




FIG. 5: Variation of Q a ct(t) with time, for different p. The 
arrows mark r a . 

sizes, until the entire space is filled up. Note that the ob- 
served largest cluster dependens on system-size and thus 
the tail of P(s) is affected by finite-size effects. 

Finally we investigate to what extent the dynamics is 
encoded in the structure. We have seen that structural 
relaxation is correlated with the development of the mo- 
bile regions. Therefore it is interesting to check how the 
properties of these regions depend on the trajectories that 
start from the same initial configuration. For this, we 
check along several such trajectories how different are the 
configurations of active sites formed at the same observa- 
tion time. To quantify that, we define the overlap func- 
tion Q act (t) = {(q aP {t)) lc - n 2 act (t))/(n act (t) - n 2 act (t)), 
where q a/3 (t) = £\ n?(i)n<?(t), with m(t) = 1 if site 
i is active at time t and rii(t) — otherwise. L)i C is 
the average over the isoconfigurational ensemble [3J, i.e. 
the ensemble of all possible trajectories starting from the 
same configuration of which a and (3 are two different 
members. Defined in this way, Q ac t(t) = 1 if at time t 
the configuration of active sites for two different trajec- 
tories are exactly the same and Q ac t(t) = if the two 
configurations are totally different, apart from the trivial 
statistical overlap. Thus Q ac t is a direct measure for the 
influence of the initial structure on the mobility regions. 

In Fig. [51 we show the time dependence of Q ac t{t) 
for different p. At short times Q ac t(t) is independent 
of p and its value is small since the mobile particles can 
find random directions to explore, which result in differ- 
ent configurations of active sites and hence a small Q ac t- 
Subsequently Q ac t(t) has a peak at a time which approx- 
imately corresponds to the time at which n act (t) enters 
the final regime of growth, see Fig. [3l and which is much 
smaller than r Q (marked by arrows). For large value of 
p the height of this peak is close to unity and it becomes 
extremely broad so that Q ac t(t) is quite large even at r a . 
This shows that the configuration of active sites is, in a 
broad time window, basically independent of the trajec- 
tory and thus encoded in the initial structure, even for 
times of the order of r a (but depends of course on the 
initial configuration) . Note that these results seem to be 
in contradiction to the claims made in Ref. Q since there 
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FIG. 6: The time evolution of overlap, As(t), for different p. 
The arrows mark r Q . Inset: Comparison of different temporal 
functions at p — 0.88. 

it was argued that the structure influences the dynamics 
only on time scales much shorter than r a . However, in 
that work the authors considered this influence on the 
level of individual particles whereas the overlap Q ac t(t) 
considered here is a collective quantity. Therefore there 
is not necessarily a contradiction. 

What structural property determines mobility? It is 
obvious that for a site to be active, it needs empty space 
in its neighborhood. Earlier work on lattice models has 
shown that the high propensity sites are located near 
clusters of empty sites [13] • In order to investigate the 
relation between structure and mobility we define an- 
other overlap function. A lattice site is considered to be 
a generalized vacancy of type k if the site and its six near- 
est neighbors contains a total of at least k holes. Then 
we calculate, as a function of time, the overlap of the 
active sites with these generalized vacancies. The quan- 
tity we measure is A k (t) = £\ afni(t)/J2i n i{t) ~ « fe - 
Here, a!? = 1 if at t = the site i is a generalized va- 
cancy of type k and otherwise. Thus, A/-(f) is the 
probability that a generalized vacancy of type k in the 
configuration at t = is an active site at time t, with 
the trivial overlap a k (the density of generalized vacan- 
cies of type k in the configuration at t = 0) subtracted. 
Such a quantity enables us to have a good measure of 
the correlation between structure and dynamics. For low 
density, i.e. p = 0.65, A${t) has a fast decay, see Fig. |6] 
(Other value of k have a qualitatively similar behavior.) 
With increasing p, the characteristic time-scale for the 
decay increases with the tail becoming stretched in shape, 
the underlying slowing down of the relaxation process re- 
sulting in the retention of the memory for longer time. 
In the inset of Fig. [HI we have plotted different correla- 
tion functions at p = 0.88: The overlap function As(t), 



the self- intermediate scattering function F s (ir, t) and the 
fraction of active sites n act (t). At this density, the a- 
relaxation time r Q w 4.4 x 10 8 . For this time-scale, the 
fraction of active sites is n act (T a ) ~ 0.50. However, if 
one extrapolates A 3 (t) to these time-scales, the overlap 
is small. Therefore, although measurement of Q ac t(t) 
showed that the configuration of active sites, at t ~ r Q , 
is significantly determined by initial structure, the use 
of generalized vacancies does not fully demonstrate this 
strong dependence. Hence, a better characterization of 
initial structure is necessary for improving prediction for 
formation of mobile regions. 

In conclusion, we have shown that, for the KA lattice 
gas, the a— relaxation occurs via the percolation of mo- 
bile regions in which particles move cooperatively. These 
regions are encoded in the initial structure of the sys- 
tem with a memory time that is on the order of t q . We 
emphasize, however, that the initial structure does not 
necessarily determine the trajectory of an individual par- 
ticle but only the location and the shape of the regions 
in which cooperative dynamics is observed. 
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